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Estimating the parameters that dictate the dynamics of a quantum system is an important task for quantum 
information processing and quantum metrology, as well as fundamental physics. In this paper we develop a 
method for parameter estimation for Markovian open quantum systems using a temporal record of measurements 
on the system. The method is based on system realization theory and is a generalization of our previous work 
on identification of Hamiltonian parameters [Phys. Rev. Lett. 113, 080401 (2014)]. 


I. INTRODUCTION 

Recent years have witnessed the rapid progress of quan¬ 
tum information processing technologies, including convinc¬ 
ing demonstrations of a quantum advantage in several appli¬ 
cations including communication and sensing. Such technolo¬ 
gies require the precise fabrication and manipulation of quan¬ 
tum degrees of freedom, and as a result, much effort is in¬ 
vested into understanding and precisely identifying the quan¬ 
tum dynamics. 

In Ref. [ 1 ] we developed a technique for identifying a pa¬ 
rameterized Hamiltonian from time traces of expectation val¬ 
ues of a small set of observables. This technique was recently 
experimentally demonstrated and validated in Ref. [2]. In the 
current paper we generalize this work to enable identification 
of parameterized open system evolution that can be described 
by a Lindblad master equation. This expands the applicabil¬ 
ity of this type of system identification approach that utilizes 
time traces of observables. As in Ref. [1] we consider only 
finite dimensional systems, and assume that the system can be 
reliably prepared in a small number of initial states, and pos¬ 
sesses observables whose expectation value can be sampled 
over a period of time. Nuclear magnetic resonance [2, 3] and 
ensembles of neutral atoms [4] are two typical examples of 
physical systems for which these assumptions are valid. 

Our approach can be considered a generalization of tradi¬ 
tional spectroscopic methods such as Ramsey interferometry 
in which spectral features of time-dependent data are used to 
infer values of underlying system parameters [3, 5]. This in¬ 
ference is simple in the case of Ramsey or Rabi measurements 
where the relation between spectral features and the parame¬ 
ters is straightforward. In more complex situations, the re¬ 
lationship can be too complex to know a priori. Moreover, 
it may not even be known whether the measurements per¬ 
formed can identify a parameter of interest. Such complex 
situations beyond conventional spectroscopy can occur even 
in small systems such as a few atoms or spins. We develop a 
systematic way to perform parameter estimation in such com¬ 
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plex situations, providing a criterion of whether the unknown 
parameters are estimable given the set of measurements avail¬ 
able, and if so, prescribing a data-driven algorithm to identify 
them. 

In addition to Ref. [ 1 ], previous studies that have examined 
this kind of generalized spectroscopy to estimate Hamiltonian 
or open system parameters from time traces include Refs. [6- 
18]. We draw particular attention to Ref. [15], which studied 
the problem of identifying open quantum systems in the same 
setting that we do. The solution constructed in that work is 
related to the one that we will present below in the sense that 
both consider observable time traces in the Laplace domain 
and attempt to solve equations encoding the relation between 
the unknown parameters and the measured signal. The prin¬ 
cipal difference is that whereas Ref. [15] expresses the un¬ 
known parameters directly in terms of the measured signal, 
we use model realization theory to first construct a minimal 
model of the system from the measurements, and then esti¬ 
mate the parameters from this minimal model. 

The remainder of this paper is organized as follows: in sec¬ 
tion II we formulate the system identification problem in the 
Markovian open quantum systems context. Then in section III 
we describe the identification algorithm based on model real¬ 
ization, and highlight some of the distinctions between this 
algorithm and the corresponding algorithm for closed systems 
developed in Ref [ 1 ]. We illustrate the open system algorithm 
with an example in section IV, and in addition, we present sev¬ 
eral case studies in Appendix A that exemplify the types of 
equations that must be solved in order to identify the parame¬ 
ters. Finally, section VI concludes the paper with a discussion 
of results and future directions. 


II. PROBLEM FORMULATION 

The problem of interest is the identification of unknown 
parameters that dictate the dynamical behavior of a quantum 
system. The system state is represented by a density matrix 
p £ where p'^ = p, p > 0 (p is a positive semidefi- 

nite matrix), and tr p = 1. We assume the system dynamics is 
governed by a quantum Markovian master equation in Lind- 
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bladform[19, 20]: 

p = Cp= ChP + CdP, (1) 

where Ch is the generator of Hamiltonian dynamics 

CHP = -i[H{e),p] (2) 

with 9 a vector of unknown parameters (h = 1), and C£> rep¬ 
resents a general Lindblad dissipative generator 

= 2 5,^ [[Fj,pFl] + [F,p, Fl]) . (3) 

Here {iFk} is an orthonormal basis for the Lie algebra 
su{N), and the Hilbert-Schmidt inner product is defined as 

{iFm^iFn) = tr(F^F„). Letting G = igjk)^j,Zl^ we have 
that Gi = G > 0 [20]. Expanding the Hamiltonian in this 
orthonormal basis, we obtain a parameterized Hamiltonian in 
the form: 

E m 

h^{e)Fm, ( 4 ) 

m—1 

where hm — tr{HFm) G K are some known functions of 9. 
We take hm as the unknown parameters from now on since 
solving 9 from hm is an algebraic problem. The unknown 
parameters that we want to estimate are the Hamiltonian pa¬ 
rameters hm and the Lindblad coefficients pj^. 

First we define the structure constants Cjki € C of the Lie 
algebra su{N) with respect to the orthonormal basis 
through the commutator of basis elements: 

[iFj,iFk]=J2^^^ CjmFi. ( 5 ) 

For future use we also define the constants Dj^i € C through 
the anticommutator of basis elements [20]: 

{iFj,iFk} = FjkiiFi. ( 6 ) 

The dynamics of the expectation value of an observable F^, 
i.e., Xn{t) = tr Fnp{t), can be written as 

in{t) =tlFn{CHp{t)+CDp{t)). (7) 

From [1, 19], we derive that 

~ {Qnp F Rnp) Xp{t) F (^) 

where 

E m 

, Fmnpb'mj 
m—1 

Rnp = — T gjk{Cnjl{Cklp + Dklp) 

^ .j. (y) 

F Cknl{Cljp F DijpZj^ 

1 

^ R^{gik)Cnjk- 


Collecting the x„ in a vector x € ^, we obtain a linear 

differential equation describing the complete dynamics: 

x(f) = Ax(f) -I- b, x„(0) = trF„p(0), (10) 

where the matrix A S -i)x(Ar -i) elements A„p = 

Qnp F Rnp 7 and the vector b G has elements bn- 

The vector x(f) is often called the coherence vector, and is 
a complete representation of the quantum state [19]. Eq. (10) 
gives an explicit description of the quantum dynamics as a 
linear time invariant (LTI) system, and hence it enables appli¬ 
cation of results from classical linear systems theory. 

Next we turn to the observables being monitored. The mea¬ 
surement output consists of a vector of time-dependent ob¬ 
servable expectation values: 

y(f) = [(Oi(f)), ( 02 (f)), (11) 

In most physical systems one only has access to a limited set 
of observables {e.g., local observables in a many-body sys¬ 
tem). However, as a result of the dynamics, it is possible that 
all the parameters defining Eq. (1) are imprinted in the time 
evolution of the monitored observables. This is the notion we 
wish to formalize and exploit. We begin by expanding each 
monitored observable Oi in the basis {fF„} of the Lie algebra 

su{N), i.e., Oi = 'Zln=i On Fn. With this expansion we can 
define the output vector as 

y(f) = cx(t) (12) 

(•i) 

with the i-th row of the matrix c being the elements On ■ Also 
define the set Af = { F^^ , F^^ i • j F^^^ }, where iz is a vector 
of length p, as the collection of unique basis elements that 
appear in the expansion of all the measured observables. 

We want to use the dynamical equation governing the time 
evolution of only the observables being monitored to estimate 
the unknown parameters. It is possible that the evolution does 
not couple the elements in Ai to all elements of the basis. This 
is equivalent to A possessing block diagonalizable structure 
and the elements of A4 being coupled only through a proper 
subblock. A constructive procedure to find the relevant basis 
elements that couple to the measured observables is the fol¬ 
lowing generalization of the filtration procedure in geometric 
control theory [21]. First we define the adjoint generator of 
dynamics through in = ti{Fn{Cp)) = tr((£'iF„)p). Explic¬ 
itly, A = F C^oX, with 

C^hX = -i[X,H{9)\ 

and 

Fix = i Y.Tk=l 9,k {[F,,FlX] F [XF„Fl]) . 

Then for the Lindblad dynamics in Eq. (1) define an iterative 
procedure as 


Go = M, 

G, = 


(13) 
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where 


A. Eigenstate realization algorithm 


= {Fj : tY{FjC'^g) ^ 0, where g e G,_i}. 


That is, at each iteration we compute the adjoint evolution of 
each of the elements of Gi-\, and if the result has nonzero 
inner product with a basis element not already in Gi_i, this 
basis element is added to G^. Since the dynamical system is 
finite dimensional, this iteration will saturate at a maximal set 
G after finite steps, and we refer to this set as the accessible 
set. Now writing all the Xk with F^ G G in a vector of 
dimension K < iV^ — 1, the dynamics for this vector is given 
by 


d X C 

— Xq = AXa -f b, 

dc 

(14) 

y{t) = CXa(f), 

(15) 


where Aisa K x K sub-matrix of A, b is a iT x 1 sub-vector 
of b, and c is a p x if sub-matrix of c; i.e., all contain only 
the elements necessary to describe the evolution of the subset 
of observable averages collected in x^, and how these define 
the measurement traces. 

Finally, we assume that the system is prepared in a fixed, 
known initial state x(0), and the corresponding initial state 
forEq. (14) is Xa(0). 


III. IDENTIFICATION ALGORITHM 


Suppose that we can measure the expectation values of 
some observables in A4 at regular time instants jAt, where 
At is the sampling period. Denote the measured expectation 
values of the observables as {y(jAf)} [31], which is the out¬ 
put of the following discretized form of Eq. (14): 


Xa(j + 1) =Ad 


Xa(j) 


/•0 + l)At 

JjAt 


dre ^’’b 


(16) 


y(j) =CXa(j), 


The first stage of the estimation algorithm is to construct a 
minimal realization of the system based on input/output infor¬ 
mation. This is achieved by ERA in three steps as follows. 

Step 1: Collect the measured data into an rp x s matrix 
(generalized Hankel matrix) as: 


y{k) y{k + ti) ■■■ y{k + ts-i) 

yUi+k) y{ji + k + ti) ■■■ y{ji+k + ts-i) 

-y{jr-l+k) yijr-l+k + ti) ■■■ yUr-l+ k + ts-l). 

with arbitrary integers ji(i = 1, • • •, r — 1) and ti (I = 1, ■ ■ ■, 
s - 1). 

Step 2: Eind the singular value decomposition (SVD) of 

Hrs(O) as 

H,,(0) = P 


= [Pi P2 


E O' 

'QT 

0 0 

Ql_ 


where P G ^rpxrp^ Q S are both orthonormal, and E is 
a diagonal matrix with the non-zero singular values of (0) 
determined up to numerical accuracy e, i.e., Tin > e for all 
i < riY. where ns is the dimension of T. The matrices Pi, 
P 2 , Qi, Q 2 are partitions with compatible dimensions. When 
the measurement time traces are noisy, determining the cutoff 
parameter e (and hence ns) can be difficult since the noise 
can lead to non-decaying singular values. In this case, one 
can choose ns by demanding that the realization produced 
by ERA is of the correct order. We will return to this issue 
when we specify how the ERA realization will be used (see 
discussion after Eq. (18)). 

Step 3: Eorm a realization of the system (16) as Ad = 
T-^ Pj'ilrsil)QiT-^ ,c = ETpiSiandx(O) = T^Qjci, 
where Ej = [Ip, Op, • • • , Op] and ei is the first column of Ig. 
The triple (A,i,c,x(0)) forms a realization of Eq. (16), and 
since the output data is an invariant of realizations, this triple 
generates the same output as the original system: 


where Ad = and for brevity of notation we use 

xa(i) = Xa(jAf) and y(j) = y(jAf). 

Eq. (16) defines a discrete time LTI system and we will now 
use invariants of different realizations of an LTI system to 
identify the unknown Hamiltonian and Lindblad parameters 
that generate the dynamics. To this end, we need to construct 
a realization from the measurement time traces. There are 
many methods for constructing a realization of a linear dy¬ 
namical system from input-output data in linear systems the¬ 
ory [22]. In [1], we presented and utilized a method called 
the eigenstate realization algorithm (ERA) [23] for identify¬ 
ing Hamiltonian dynamics, and in the following we show that 
this method can be used in this case of open system dynam¬ 
ics in Lindblad form also. Eor completeness, we include the 
specification of ERA here. 


y(j) = cA^x(O), for all j > 0. (17) 

Note that although the original system response is composed 
of a response to a non-zero initial state and a response to the 
forcing vector b {i.e., Eq. (16)), the ERA realization is com¬ 
posed of only an initial state response. This is because at the 
level of input-output realizations it is not possible to distin¬ 
guish between the response to initial states and the response 
to forcing inputs, and therefore ERA lumps all responses into 
one type. 

This completes the specification of the ERA algorithm. It 
results in a realization of the discrete time dynamical sys¬ 
tem (16) in the form of the triple {Ad, c, x(0)). 

We note that the measurements results do not have to be 
uniform in Step 1 of the algorithm; in particular, if it is known 
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that some measurement results are particularly noisy or cor¬ 
rupted, those data points can be discarded by choosing appro¬ 
priate integers ji (i = 1, ■ ■ ■, r — 1) and ti (I = 1, ■ ■ ■, s — V). 
This data filtering can reduce the estimation error caused by 
measurement noise and outliers. 


B. Estimation algorithm 

In order to estimate the parameters in the original dynam¬ 
ical system we now convert the discrete-time LTI system re¬ 
alization obtained from ERA, (A^^, c, x(0)), to a continuous 
time realization, by letting A = log A^/At, where the log¬ 
arithm in this definition is the principal branch of the natural 
logarithm. The accurate conversion from a discrete time sys¬ 
tem to a continuous time system relies on the sampling time 
At being sufficiently small to capture all continuous time dy¬ 
namics. For the Hamiltonian system identification case we 
were able to specify conditions on the sampling time based 
on the Shannon-Nyquist criteria, since in this case the output 
time traces y(f) are guaranteed to be band-limited [1]. How¬ 
ever, there is no such guarantee in the Markovian open system 
case since now the time traces are generally decaying oscilla¬ 
tions, or time-limited sigmals. Therefore, we do not currently 
have conditions on the sampling time for this conversion to ac¬ 
curately provide a continuous-time realization {i.e., to avoid 
aliasing effects), but note that one can generally estimate a 
valid sampling time from knowledge of the intrinsic frequen¬ 
cies in the system. In engineering, it is typical practice to 
sample such time-limited signals at 6 or 8 times the fastest 
frequency in the signal [24], and such a heuristic suffices in 
our setting as well, as the example in section IV and the case 
studies in Appendix A bear out. Additionally, we note that 
although the Shannon-Nyquist criteria enables a formal spec¬ 
ification of the minimum sampling time in the Hamiltonian 
estimation scenario, to apply it requires some knowledge of 
the system’s spectrum. In the absence of this knowledge, the 
situation is the same in the closed and open systems cases; 
one needs to estimate a valid sampling time and possibly also 
try multiple sampling times. 

Now, to estimate the unknown parameters, we use the fact 
that the system input-output relations are invariants of differ¬ 
ent realizations. We work in the Laplace domain in the fol¬ 
lowing in order to form algebraic equations for the unknown 
parameters. By equating the Laplace transform of the outputs 
of the original system and the ERA realization we get; 

c(sl — A)“^ |xa(0)-l—I = c(sl — A)“^x(0). (18) 


and the coefficients in this expansion can be equated to yield 
polynomial equations for the unknown parameters. Solving 
these multivariate polynomial equations leads to the identi¬ 
fication of gjn and Ujk- In the case studies encountered in 
Appendix A we are able to express the resolvents exactly by 
computing the matrix inverses symbolically, but one may need 
to resort to the expansion in s in more complicated cases. We 
discuss this issue further in Section VI. 

Note that in Step 2 of ERA dictates the maximum or¬ 
der of the polynomial on the right hand side of Eq. (18). This 
suggests that we should choose ns to be the order of the de¬ 
nominator polynomial in the left hand side of Eq. (18), which 
can be calculated from symbolic computations and obtained 
as an irreducible rational function. This choice coincides with 
the rank of H^s (0) when there is no noise in the measurement 
time traces. 

We highlight a crucial difference here between parameter 
estimation in the closed system and Markovian open sys¬ 
tem scenarios. For the former, b = 0, and if the triple 
(A, c, Xa(0)) form a minimal realization (i.e., is controllable 
and observable), then the Laplace transform on the left hand 
side of Eq. (18) is gauranteed to have a canonical form as a 
ratio of polynomials Q{s)/P{s) [22], with 


P(s) = det(sl - A), 


Q{s) = det 


I 0 
0 0 


A Xa(0)l\ 

c 0 Jy ■ 


(19) 


Having this form enabled us in Ref. [1] to avoid explicitly 
computing the inverse (si — A)~^ or its expansion in a power 
series. These are both computationally intensive to compute 
since the matrix A is a symbolic matrix containing the un¬ 
known parameters. However, in the open system case, where 
the left hand side of Eq. (18) has no known canonical form, 
one cannot avoid performing this symbolic inverse or power 
series computation. This is a critical difference in computa¬ 
tional difficulty between the closed and open system parame¬ 
ter identification problems. 

We note that by converting back into the continuous time 
domain and formulating the Laplace transform relation in 
Eq. (18), we have converted the problem of estimating param¬ 
eters to one of solving polynomial equations. This is in con¬ 
trast to directly estimating parameters using Eq. (17), which 
would involve solving transcendental equations for the un¬ 
known parameters. This simplification of the equations re¬ 
lating the unknown parameters to the measured data is one of 
the primary advantages of our approach. 

We conclude with two further comments on the above esti¬ 
mation algorithm; 


This equation relates the unknown parameters to the mea¬ 
sured data through the ERA realization. Explicitly, the right 
hand side of Eq. (18) is completely determined by the mea¬ 
sured data, and the left hand side is in terms of the Hamilto¬ 
nian parameters hm and Lindblad coefficients gjk. The resol¬ 
vent expressions on the right and left hand sides of Eq. (18), 
(si—A)~^ and (si—A)~^, can be computed symbolically, or 
alternatively the expressions can be expanded in powers of s. 


1. The initial state may have to be chosen carefully. For 
example, if x^ is zero or an eigenvector of A, it leads 
to no sensitivity in the output to any of the unknown 
parameters. Care must be taken to avoid such degener¬ 
ate cases. In fact, running the algorithm with multiple 
initial states leads to more polynomial equations of low 
order and thus can help to solve for the unknown pa¬ 
rameters more efficiently. 
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2. This system identification algorithm can result in mul¬ 
tiple estimates of the unknown parameters, all of which 
satisfy the input/output relations captured by Eq. (18). 
This is because several Markovian generators can gen¬ 
erate the same map between an input state and mea¬ 
surement time trace, and hence are equivalent from an 
input/output perspective [25]. When the algorithm re¬ 
sults in multiple parameter estimates, one has to ap¬ 
peal to prior information, or add resources such as addi¬ 
tional input states or observable time traces in order to 
uniquely specify a parameter set. 

The following section explicitly demonstrates the algorithm 
developed above and Appendix A presents additional case 
studies. 


IV. DISSIPATIVE ENERGY TRANSEER 


In this section we apply the estimation algorithm developed 
above to a physically relevant example: energy transfer be¬ 
tween two qubits at hnite temperature. The qubits could rep¬ 
resent the relevant energy levels of an atomic or molecular 
system, or spin-| systems. The Markovian master equation 
for the dynamics takes the form 


dp 


^ - p) 


dt 


k^l 

I \ '' c\ — ( k k ^ k k ^ k k 

^2^^9k ( 2 ^ + 

k^l ^ 

+ 2gt (a'lpa’l - , (20) 


where the Hamiltonian for the system is given by 

+ (21) 


This is a model for two possibly detuned qubits with co¬ 
herent energy transfer dynamics, independent dephasing, 
and incoherent excitation-relaxation. The rates of excita¬ 
tion/relaxation (g'^ jg'^') are functions of the temperature of 
the environment of the qubits [26]. There are nine pa¬ 
rameters that dictate the dynamics in this model: 0 = 
{uji,uj2,6i,ji,-f2,9t,9i,92^92)- For convenience, we de- 
hne 


^i=9t +9i , Pi=9t -9i : 
’^2=gt+92> P2 =9^-92- 


( 22 ) 


Assume that the observable being measured is zi{t) = 
O'" ^ 2 {t) = (o'f(f)), or both, since we will be inter¬ 
ested in exploring the benefits of measuring multiple observ¬ 
ables. The state vector determined by the accessible set is the 
same regardless of whether one or both of those observables 
are measured, and is explicitly specified by = [(ct^), (<7^), 

whereas the dynamics of 


Xa is determined by Eq. (14) with 


A = 

'—2i2i 

0 

0 

(5i 

-51 

0 ■ 

0 

-2v2 

0 

-(5i 


0 

0 

0 

-Vs - 7s 

UJ2 

LOi 

0 

-<5i 


—UJ2 

-Vs - 7s 

0 

Wl 

<5i 


-Wi 

0 

-Vs - 7s 

W2 

0 

0 

0 

-Wi 

—0J2 

-Vs- 7s. 


b = [pi, P2, 0, 0, 0, 0]^, 


where Vs = vi + V 2 and 7s = 71 -I- 72. The matrix dehning 
the output depends on which observable is being monitored 
and is given by 

{al): = [1,0,0,0, 0,0] , or 

(a,"): c,, = [0,1,0,0, 0,0] , 

or the concatenation of these two row vectors into a 2 x 6 
matrix, if both observables are being monitored. 

Letting the initial state be |'!/)i(0)) = |0) (g) (|0) -I- |l))/-\/2, 
we have that 

x,(0) = [1, 0, 0, 0, 0, 0]^. (24) 

Since the matrix A is small, we can symbolically calculate the 
resolvents in Eq. (18), to obtain the Laplace transform of the 
two possible measurement traces zi{t) or Z 2 {t): 

Zi{s) = £21 (si - A)“^Xa(0) -h C2j(sl - A)“^b/s 
_ + 923^ + gis-I-go 

S® -I- P4S'^ + Pss^ + P2S^ +Pis 

with 

93 = Pi + 2z2i + 4j 22 -I- 27s, 

92 = + %viV2 + 2(217s -I- 2pii'i -\- 51^2 + 6i/27s 

+Ap4iy2 + 7? + 2pi7s + Wrf, 

gi = 2 diPi 26‘^P2 261121 261122 + 2i5^7s T Pi^i 

+6piV2 + Avivl + 2vfv2 + pi7s + 2:/27 s + 

6~Pi^d 2,12210^1 -\- 2 v 2 -\- 6pii2ii22 + 2pii2i^g 
+ Gpil22"fs + 4,12ll22"fs 

go = 2 piV 2 + 26 lpi^s+ 26 lp 2 'ys + 4 :Pit^ivl+ 2 piiylv 2 
+2pii22"i1 + 4p,iJ^27s + 2piV2Uj\ + 26\pii2i 
+26lpii22 + 26fp2i2i + 26fp2i22 + 4/rii^ii/27s, (26) 

and 

P4 = 4i/i-I-4 i^ 2 + 27s, 

P3 = 46f + bvl + \Av 1 V 2 + 6:^i7s -I- 5i2| -p 6^27s 
+ 7s + 

P2 = d>6‘ll>i+^6'lv2 + ^6'1'ys +2vl + I2^vlv2+ 

-1-14^1 r-l -I- \&viV2'ys + + 2viuj'^ 2r'| 

+4j^27s + 2v2^‘1 + 2v2^‘dJ 

Pi = A6\v^ i6\v\V2 + 6i6\vi"fg 46^122 -|- 4 ( 5 ^z.' 27 s 

+412^122 + 812^122 + 8l2fl22js + 4z/iz/| -I- 8121122 JS 
+4121122^“^ + 4i2i122Uj'^, (27) 
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where ujd = wi — 0 J 2 - And, 

Z2{s) = Cz2(sl - A)~^Xa(0) + C^^{sl - A)“^b/s 

ras^ + r2S^ + ris + rg (28) 

s® + + pas^ + p2S^ + Pis ’ 

with 

?'3 = M2, 

r2 = 2(5i + 4p2J^i + 2p2J^2 + 2p27s, 

Ta = ‘2S‘^pi “t" ‘26‘^P2 ~t“ 2 ^a^i “h 2 ()a ^2 “h 2^i7s ~h 
+P2^'| + Ai27s + + 6/r2l^lJ^2 + 6p2i^l7s 

+2p2J^27s> 

To = 2p2J^i + 2(5a/ri7s + 2Slp2ls + 2p2Vivl + 4:p2vlv2 
+2p2i^i7s + 4p2J^i7s + ‘2.p2ViU]l + ‘i.SlpiVi 
+2SfpiV2 + 2dlp2Vi + 2Sfp2V2 + ^P-2^lV2ls (29) 

and Pi defined as in Eq. (27). 

It is clear from these expressions that we can identify Vi, 
Pi, which then allows identification of gf, g~. At the same 
time, note that only the linear combinations 7 s = 71 + 72 
and uJd = uji — UJ 2 occur in the above equations but not the 
individual parameters 7 ^ and uJi {i = 1, 2). This implies that 
only these linear combinations can be determined from the 
measurements and the initial state, but not the individual pa¬ 
rameters that enter them. These linear combinations describe 
the energy difference between the qubits and the dephasing- 
induced broadening of this energy difference. The physical 
observables encoding average population of the excited state 
of either qubit only allow determination of these collective (in 
the case of 7 s) or relative (in the case of Wd) properties of the 
system. Furthermore, another restriction that we can imme¬ 
diately observe is that only even powers of and ojd occur 
in all of the polynomials above, and therefore we only expect 
to determine these parameters up to a sign difference. The 
signs of these parameters are not estimable by the local mea¬ 
surements we have chosen, and prior knowledge, or additional 
initial states or measurements are necessary to identify these 
signs. 

Before we present the results of the estimation simulation, 
we outline four possible modes of estimation that could be 
performed given that there are two possible observables: 

1. In mode 1 , only the measurement trace Zi(f) is used 
to construct the Hankel matrix and system realization 
(c = c^j^), and the polynomial system used to solve 
for the unknown parameters is constructed using only 

Zi{s). 

2. In mode 2, only the measurement trace Z 2 {t) is used 
to construct the Hankel matrix and system realization 
(c = Czj), and the polynomial system used to solve 
for the unknown parameters is constructed using only 

Z 2 {S). 

3. In mode 3, the measurement traces zi {t) and Z 2 (t) are 
used to construct the Hankel matrix and system realiza¬ 
tion (c = Czjzj)^ ^nd the polynomial system used to 



FIG. 1: (Color online) Measurement time traces of (^a\{t)) and 
(erf (t)) for the energy transfer example presented in section IV. The 
initial state and system parameters are specified in the main text. 

solve for the unknown parameters is constructed using 
only Zi{s). 

4. In mode 4, the measurement traces zi{t) and Z 2 {t) are 
used to construct the Hankel matrix and system realiza¬ 
tion (c = 62 ^ 22 )^ ^nd the polynomial system used to 
solve for the unknown parameters is constructed using 
only ^ 2 ( 3 ). 

Obviously, modes 1 and 2 are practically more attractive 
since they only involve collecting one observable’s time trace. 
However, there may be a benefit to consider modes 3 and 4 
since in these modes of estimation one is providing the algo¬ 
rithm with more data with which to construct the realization 
(although the accessible vector does not change, and so some 
of this additional data is redundant). One could also imagine 
a fifth mode where one uses the measurement traces zi (t) and 
Z 2 {t) to construct the Hankel matrix and system realization 
(c = then constructs the polynomial system from 

the definitions of both Zi{s) and Z 2 {s), i.e., some polyno¬ 
mial equations from the definition of the coefficients in Zi (s) 
and some from the definition of the coefficients in ^ 2 ( 3 ). For 
our example, where the number of parameters is small enough 
such that one can get enough equations from using the defini¬ 
tion of just one of the Laplace transforms, we did not find any 
advantage to using this fifth mode of estimation, and so do not 
investigate it further. 

To illustrate the algorithm we fix the nominal parameters of 
the system as wi = 1.3MHz, uj 2 = 2.4MHz, 5i — 0.5MHz, 
7 i = 0.03 MHz, 72 = 0.035MHz, g+ = 0.02n(wi)MHz, 
g- = 0.02(fi(wi) -f l)MHz, with 

n{uj) = ^ ^ - 

gfcijT _ I 

being the Bose-Einstein distribution at temperature T. In the 
following we fix the temperature at fegT = O.Swi. For a sys¬ 
tem with these parameters, we generate a time trace of zi {t) 
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and Z 2 {t) from t = Q lo t = ti = 60/is, with Ai = O.Ol^s. 
The resulting time traces are shown in Fig. 1. We form the 
Hankel matrix using either one or both of these time traces, 
depending on the mode of estimation, with ti = ji = i 
(i.e., using every data point), and r = s = ^^1^^ . The 
ERA realization of the discrete time system is formed with 
ns = rank(Hj.s(0)), and the corresponding continuous time 
system is formed using the prescription given in section III B. 
Finally, the Laplace transform expression on the right-hand- 
side of Eq. (18) is computed to obtain [32] 

s'* -f 0.2702 s 3 -f 1.7302 s 2 -f 0.072s - 0.0034 
* ~ s5 + 0.3624?’ + 2.2569s3 + 0.3243s2 -f 0.011s 

and 

-0.0176s 3 -f 0.4944 s 2 -f 0.0209s - 0.0039 
^ ~ s5 -f 0.3624s4 -f 2.2569s3 -f 0.3243s2 -f 0.011s' 

Combining these expressions with corresponding Laplace 
transforms Eqs. (25) and (28), we obtain a system of poly¬ 
nomial equations for the unknown parameters (we can choose 
the seven simplest equations since there are seven unknown 
parameters). We solved these equations using PHClab, the 
Matlab interface to the PHCPack libraries for solving poly¬ 
nomial systems [27]. 

Table I shows the results of solving for the parameters when 
the measurements are noiseless. We found no difference be¬ 
tween modes 1 and 3, and also between modes 2 and 4, 
when the measurements are noiseless, and therefore we only 
present results from modes 1 and 2 of estimation in table 
I. For mode 1, the first observation is that there exist sev¬ 
eral estimates that deviate from the nominal values as shown 
on the first line. This is because in this case we choose the 
lowest order polynomials defined by the coefficients of Zi{s) 
and the resulting polynomial system has multiple solutions, 
among which the nominal set is only one of them. Secondly, 
as expected from the above observation that only even orders 
of i5i and ojd enter the polynomial system, the signs of these 
parameters are indeterminate. In estimation mode 2, the es¬ 
timate quality increases and in fact, the only uncertainty is in 
the sign of 5i and oJd- This is an example of how the choice 
of observable dictates the quality of the estimation. 

In summary, we see that given noiseless measurement 
records, the above estimation algorithm can determine the un¬ 
known Hmailtonian and Lindblad parameters to within the 
limitations of the data (e.g., in the above example, the limi¬ 
tations were that uii, uj 2 , 7 i , and 72 are not individually es¬ 
timable and that the signs of and uid are not estimable). 

V. NOISY MEASUREMENTS 

In this section we investigate the robustness of the sys¬ 
tem identification algorithm by reexamining the two-qubit 
dissipative energy transfer example with noisy measurement 
traces. In principle, the noise on expectation values of ob¬ 
servables can be made arbitrary small since the signal-to-noise 
decreases as 1 /y/N, where N is the number of measurements 


that are averaged to estimate the expectation value of the ob¬ 
servable. Because of this, measurement noise is especially 
small in ensemble systems like NMR [2]; however, in sys¬ 
tems without natural access to ensembles, e.g., a single super¬ 
conducting qubit, noise on time traces cannot be neglected in 
practical situations and we must assess the robustness of the 
above system identification algorithm to noise. 

We consider the same system as in section IV, with ac¬ 
cess to time traces of one or both observables: and 

^cr^(f)). Suppose that these time traces are corrupted by ad¬ 
ditive Gaussian noise 

where ^i(j) ~ N{0,(J^) for all i,j. Since the expectation 
values are formed by averaging many independent measure¬ 
ment outcomes, this is a reasonable model for the noise by the 
Central Limit Theorem. We construct the Hankel matrix and 
perform the estimation in exactly the same way as section IV, 
with the only difference being that ns is fixed to be 5 instead 
of the rank of the Hankel matrix Hrs(O), since the order of 
the denominator polynomial on the left hand side of Eq. (18) 
(which takes the form in Eqs. (25) or (28) for this example) is 
5. 

To assess the quality of estimation, we compute the relative 
error in estimation for each of the seven parameters as 


where Oi and 9i are the nominal and estimated values of the 
parameter, respectively. If the estimation produces multiple 
solutions, then we choose the one with the least sum of er¬ 
rors X]I=i 6*- We generate M = 500 instances of noise with 
given standard deviation a, and calculate the estimation errors 
for each instance. These errors are then averaged to yield a 
mean relative error Cj, which captures the performance of the 
algorithm. We evaluate these mean relative errors for standard 
deviations a = 0.05, 0.10, 0.15, as well as for the noiseless 
case (cr = 0). The initial state in all instances is the same as in 
section IV, and the only difference in parameters is that the du¬ 
ration of the measurement time traces is longer: = 120/j,s. 

The first observation from this simulation (data not shown) 
is that there is a difference between modes 1 and 3 of estima¬ 
tion (and modes 2 and 4) when the measurements are noisy; 
it is beneficial to use data from both measurement time traces, 
^noisy(^) form the realization even if the poly¬ 

nomial system is formed from the Laplace transform of one 
of the time traces. Therefore we present only the estimation 
modes 3 and 4. 

Figure 2 shows the mean relative error for each parame¬ 
ter as a function of the standard deviation of the measure¬ 
ment noise. The left (right) column plots under mode 3 
(mode 4) of estimation. This figure shows that the aver¬ 
age error in estimation is small for small tr, but quickly be¬ 
comes quite large. An interesting feature is that the perfor¬ 
mance can be very different under mode 3 and mode 4 of 
estimation, with one performing better for some parameters 
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Parameters 

LOd 

<5i 


V2 

Ml 

M2 

7s 

Nominal values 

-1.1 

0.5 

0.0361 

0.022 

-0.02 

-0.0176 

0.065 

Mode 1 

±1.0973 

±0.5029 

0.0677 

-0.0096 

0.0432 

-0.0815 

0.065 

±1.1 

±0.5 

0.0361 

0.022 

-0.02 

-0.0176 

0.065 

Mode 2 

±1.1 

±0.5 

0.0361 

0.022 

-0.02 

-0.0176 

0.065 


TABLE I: Estimates derived from noiseless measurement traces. The modes of estimation (mode 1 and mode 2) are explained in the main 
text. 



a a 


EIG. 2: (Color online) Average relative error of parameter estimates 
as a function of the standard deviation (ct) of the additive noise on 
the measurement traces. In all the plots, the y-axis shows the relative 
error as a percentage and the x-axis is u. The left (right) column 
shows the error for estimation under mode 3 (mode 4). The red 
circles denote the result of the better estimation mode for that param¬ 
eter (lower maximum error), and the error bars indicate the standard 
error for the average (which is taken over M — 500 instances of 
noise). 


and worse for others. Also, the performance can vary sig¬ 
nificantly among parameters. Estimates of the Hamiltonian 
parameters, uJd and ^i, exhibit the greatest robustness to mea¬ 
surement noise, while estimates of the open system parame¬ 
ters are more sensitive. In particular, estimates of pi,p 2 and 
V 2 , which are the parameters of smallest magnitude, suffer the 
most from measurement noise. It is clear that this open system 
identification algorithm is not as robust as the Hamiltonian 
version [1]. However, in the regions of low noise (u < 0.1) 


relevant for NMR and other ensemble experimental platforms, 
the performance is acceptable, especially for the Hamiltonian 
parameters, which can be estimated well even in the presence 
of dissipation and dephasing. 


VI. DISCUSSION 

We have extended the quantum system identification ap¬ 
proach developed in Ref. [1] for Hamiltonian systems, to 
Markovian open quantum systems. The approach proceeds 
by forming a realization of the quantum system from (time- 
dependent) input-output data and then uses this realization to 
form a system of polynomial equations for the parameters that 
define the system. The strengths of the approach are its abil¬ 
ity to incorporate prior information and its ability to produce 
parameter estimates even with time dependent measurements 
of only a few observables. 

As with the Hamiltonian case dealt with in Ref. [1 ], having 
access to time-dependent data enables us to directly estimate 
the generator of dynamics, which has a significant advantage 
in that it is specified by fewer parameters than the dynamical 
map at a fixed time {e.g., Kraus map or unitary). In contrast to 
the Hamiltonian case, the Markovian open system identifica¬ 
tion scenario ideally requires the symbolic computation of a 
resolvent as expressed in Eq. (18), which can be computation¬ 
ally expensive for large systems. Alternatively, as mentioned 
in Section IIIB one could avoid the resolvent computation 
by using a power series expansion of the resolvent. This ap¬ 
proach is equivalent to using the Markov parameters of an LTI 
system as the model realization invariant, as opposed to the 
Laplace transform of the output. Although this is computa¬ 
tionally advantageous, we have found that the Markov param¬ 
eters are more susceptible to measurement noise and therefore 
this approach is expected to yield less robust estimates of the 
parameters. As a consequence, we expect that symbolic com¬ 
putation of the resolvent will be the fundamental limitation 
to applying this approach to identification of large Markovian 
open quantum systems. However, as demonstrated in Sections 
A and IV, the required symbolic computation can be easily 
performed for small open systems, which are still difficult to 
identify using other approaches. 

The noise robustness results in V suggest that a direction 
for future work is to understand exactly why this realization- 
based system identification algorithm is more robust to mea¬ 
surement noise in the Hamiltonian (closed system) case than 
in the open system case. One clue as to why this may be is 
that Markovian open system parameters almost always dic¬ 
tate exponential rates of decay of measurement traces, and 
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parameter fitting to decaying exponentials is a notoriously 
ill-conditioned problem [28]. More work is required to un¬ 
derstand the properties of this algorithm under measurement 
noise and to increase its robustness in the open system sce¬ 
nario. 
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Appendix A: Case studies 


In this appendix we present some examples that illustrate the Markovian open system identification algorithm developed in this 
paper. The system studied in all examples is a ID chain of n-qubits, with closed system dynamics governed by the Hamiltonian 

(Ai) 

k=l k=l 


This Hamiltonian is often used as a model for a spin “wire” that enables quantum state transfer [29], and is a relevant to resonant 
energy transfer mechanisms between atomic or solid-state systems, e.g., [30]. 

In Ref. [ 1 ] we analyzed the identihcation of the Hamiltonian parameters of this closed system when the observable being 
directly measured is (cr^). When the generalized Pauli operators is chosen as the basis of operators {Fk}, the accessible set 
was found to be G = U {2“”/^tT^ • • • Then denoting the average 

of the observables in G as Xq = [xi,yi, ...,Xn,ynV, where xi = (crl) ,yi = (crl) and Xk = {al ■ ■ ■ ,yk = 

fork > 2, the dynamics of is determined by 


where the system matrix A is given by [ 1 ]: 


Xa 


= Ax 


aj 


A = 


0 0 —i5i 

—oji 0 i5i 0 0 

0 —0 W2 0 

(5i 0 — UJ2 0 

0 ■■. ■■■ 

0 0 

Sn-1 


0 


— Sn-1 

0 ^n -1 0 

^n—1 0 

0 -Wn 0 


(A2) 


(A3) 


Einally, in this basis the direct observation of (cr^) corresponds to c = [1, 0,0,..., 0]. 

In the following we will consider various types of open system dynamics for this ID spin chain and apply the algorithm 
developed above to estimate the parameters defining the Hamiltonian and the open system dynamics. 


1. Independent dephasing 

We first consider a common decoherence model for spin-^ systems, namely, independent Markovian dephasing on each 
spin [26]. This is described by the Lindblad terms 

n 

^ 27 fe(cr^pCT^-p), 


(A4) 
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where the superscript k indexes the spin, and the factor 2 is for the ease of expression. Comparing Eqs. (3) and (A4), we see that 
the independent dephasing model is equivalent to the choice 


gik 


Ik^ik, if tr(Xfecr^) ^ 0; 

0, otherwise, 


(A5) 


meaning that the Lindblad coefficient matrix G is diagonal and real. Again, assuming that the directly measured quantity is 
(cr^), the system matrix A now becomes: 



-71 

CUl 

0 

-h 

-Wi 

_7i 

5i 

0 

0 

-5i 

-72 

U)2 

5i 

0 

—UJ2 

-72 


0 ■■■ ■■■ 

0 


0 

0 

0 

— Sn-l 

• 1 ^n—1 0 

0 ^n—1 'Jn 

^n—1 0 Tn 


(A6) 


Comparing Eqs. (A3) and (A6), we observe that the dephasing terms introduce nonzero diagonal elements in A, whereas the 
accessible set G remains the same. In addition, for these dynamics, b = 0. All parameters in the Hamiltonian Eq. (Al) and 
Lindblad coefficients in Eq. (A4) appear in A, and therefore it is possible to identify these parameters by measuring {al). 
Consider the case of three-qubit (n = 3), whereby 



■ -71 

UJl 

0 

-5i 

0 

0 


-UJi 

-71 

(5i 

0 

0 

0 

A = 

0 

(5i 

-<5i 

0 

-72 

— UJ2 

022 

-72 

0 

52 

-52 

0 


0 

0 

0 

-52 

-73 

W3 


0 

0 

<52 

0 

-W3 

-73 


Now, if the initial state is given by |'0(O)) = (|0) + \\))/y/2 ® |00), then we have 

xjO) = [1, 0, 0, 0, 0, 0]^. (A8) 

This is an example where the Laplace transform of the measurement trace xi {t) takes the canonical form given by Eq. (19), and 
we obtain 

^ + q2S^ + qiS + gp 

s® + + P2S^ + pis + po ’ 

where the expressions for pk and qk in terms of Hamiltonian parameters and Lindblad coefficients are : 


Pa = 2(71 + 72 + 73), 

P4 = 26 l + 2S2 + 7i + 47172 + 72 + 4(71 + 72)73 +73 +^i +1^2 + <^31 

P3 = 2(7472 + 717I + 7 i 73 + 4717273 + 7273 + 7173 + 7273 + <52(271 + 72 + 73) + <5?(7 i + 72 + 273) 

+ 72W? + 73W? + 71^2 + 73^1 + (71 + 12)ujI), 

P2= 6^ + 62+ 7 i 72 + 47 i 7273 + 4717I73 + 7 i 73 + 47 i 7273 + 7273 + 72^1 + 47273W? + 73 <^i + 7 i ‘^2 + 47173^1 
+ 73<^2 + ^1^2 + (71 + 47172 + 72 + <<^1 + ^2)^3 + 2^2(71 + 7273 + 271(72 + 73) + LOi — 022^3) 

+ 2(54 ((^2 + 7 i 72 + 2(71 + 72)73 + 73 - ^1^2 + wl), 

Pi = 2(5|7 i -I- <5473 + 7273(717273 + 71(72 + 73) + (72 + 73)^?) + 73(71(71 + 73) + w?)w| 

+ (72(71(71 + 72) + W?) + + ^i( 5|(7 i + 73) + 73(27172 + (71 + 72)73 - 2a;icu2) 

+ (71 + 72 )^ 3 ) + ^ 2 ( 71(72 + 73 ) + (72 + 73)w? + 271(7273 - UJ 2 UJ 3 ))), 

Po = (< 5 l + 2(54(7172 — OJ1UJ2) + (71 + Wi )(72 + <^2))(73 + <^3) + 252((54 (7173 + UJ1LO3) + ^2(71 + <*^i) 

+ ( 7 i + Wi )(7273 - W2W3)), 
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and 

94 = 7 i + 272 + 273, 

93 = + 26 l + 7I + 47273 + 73 + 271(72 + 73) + + ujI, 

92 = 2 ( 5 |( 7 i + 72 + 73) + ^?(72 + 273) + 273(72(72 + 73) + w|) + 272 w| + 71(7! + 47273 + 73 + ^2 + ^3); 

9 i = ^2 + 7273(7273 + 271(72 + 73)) + 73(271 + 73)^1 + (27172 + 72 + w|)a;| + 2^^(7273 + 71(72 + 73) - W2W3) 
+ ^1(^2 + 27273 + 73 + ^ 1 )^ 

9 o = ^271 + <^2(^173 + 2717273 - 27iW2a;3) + (^^72 + 71(72 + + Wg). 

The Hamiltonian and Lindblad parameters can be obtained by solving these polynomial equations. 


2. Independent relaxation 

Next we consider another common form of Markovian decoherence, i.e., independent relaxation [26]. This amounts to the 
following Lindblad terms 

- \aXatp-\palat), (A9) 

k 


where is the relaxation rate of qubit k. 

Consider for example, the two-qubit case where we directly observe (cr^). Performing the generalized filtration procedure 
prescribed by Eq. (13), the accessible portion of the coherence vector is found to be Xa = [{o']:), {<^y), (c^). ('^y)’ 

{(TyCj].), (cr^cr^), (CT^ciy)]^. Comparing with Eq. (A2), we see that the Lindblad relaxation term has introduced four new elements 
into the accessible state vector, namely, (ct^), (c^), (tr^CT^), and (ct^ct^). The dynamical equation is described by Eq. (14), with 


-9i 

OJi 

0 

0 

0 

0 

0 

-h 

-UJi 

-9i 

0 

0 

0 

0 

<5i 

0 

0 

0 

-92 

a;2 

0 


0 

0 

0 

0 

— 0 J 2 

-92 


0 

0 

0 

~292 

0 

0 

-5i 

-292 - 9i 

Wl 

0 

0 

0 

~292 


0 

-Wi 

- 29 ^ - 9i 

0 

0 

0 


-29i 

0 

0 

0 

-29r - 92 

W2 


0 

0 

-29r 

0 

0 

- 0 J 2 

-29r - 92 . 


and b = 0. Note that the constant forcing vector b for the dynamics of the accessible set can be zero even if b 7 ^ 0 (as in this 
case). 

Suppose we prepare the initial state |'0(O)) = (|0) + |l))/-\/2 ® |0), then 

x,( 0 ) = [ 1 , 0 , 0 , 0 , 1 , 0 , 0 , 0 ]^. 


The resolvent on the left hand side of Eq. (18) can be computed symbolically in this case, and consequently the Laplace transform 
of the measurement trace xi is 


Xi(s) 


+ 96'S^ + ■ ■ ■ + 9is + 9o 

s® + + Pes® + • • • + pis + po 


(All) 


The expressions for the coefficients in Eq. (All) are quite involved and so we only present the most concise ones here; 


96 = TPi + 8^2 , 

95 = 26(52 + 4892 9i + 3^1 + 2tLi| + 19(5]^ + wj, 

54 = 1309]"( 5^)2 + 175^9^ +5a;?9f ^lOujlg^ + 12a;|9^ + 108(9f)^9^ + 44 ( 9 ^)®-f 185^9f + 25(9f)® + 4w^9^, 

and 


Pi = 891 + 892 , 

Pe = 2^2 + 26 ( 9 ]^ + 26(92 + 5692 9 i + SwJ + ASf. 


(A12) 





12 


These five equations and the realization formed from the measurement data can be used to solve for the five unknown parameters 
in this model. 

To treat a case where b ^ 0, we now alter the setup to consider direct measurement of a time trace of zi = (cr^). In this 
case, with the Hamiltonian given in Eq. (Al) with n = 2 and independent relaxation as prescribed in Eq. (A9), the accessible 
set becomes G = {cr], cr^, cr^a^, a^ay, cr^CT^, The dynamics of Xa is determined by Eq. (14) with 


-2gi 

0 

0 

<5i 

-<5i 

0 

0 

— ‘^92 

0 

-^1 


0 

0 

0 

-97 

U)2 

Wi 

0 



—0J2 

-97 

0 

UJi 


-<5i 

-Wl 

0 

-97 

0J2 

0 

0 

0 

-Wi 

—UJ2 

-97 


b = [-g^ , -p 2 , 0 , 0 , 0 , O]^ , 

where gj = gi + In this basis c = [ 1 ,0,0,0, 0, 0]. Let the initial state be |V'(0)) = (|0) + |l ))/-\/2 0 |0) as before, in 
which case 

xa(0) = [0, 1, 0, 0, 0, 0]^. 

Again, we can symbolically calculate the resolvents in this case, and obtain the Laplace transform of the measurement trace zi 
as 


with 


and 


^i(s) = c(sl — A) ^Xo(O) + c(sl — A) ^b/s 
_ qas^ + 925^ + qis + go 

S® + PiS'^ + PsS^ + + Pis ’ 


(A13) 


93 = -gf) 

g2 = 25^-2(5f)2-4pj-5^, 

91 = -9i {{gT+ ^9i92 + H92 f + ^d) , 

go = -2(5?(5r +92?- ‘^9i92{{9i +92? +^d)^ 


Pi = 4(gi +52 ): 

Ps = 4<5i + 5(5;^ + lAg^ g^ + 5(^2 ? + 

P2 = 2(5;^ + 52 )(4(5i + (5i ? + 651 52 + (92 ? + ^d)’ 
Pi = 4i5?(5r +9^)^ + 45r9^((9r +92? +^^1), 


(A14) 


(A15) 


where tOd = uji — UJ 2 - An interesting aspect of this example is that from Eqs. (A14) and (A15) we can identify g^, 5 ^, and 61 , 
but only uji — uj 2 - The individual transition energies of the qubits do not influence the measurement trace, and only the their 
difference does. ._ 
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